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1.  Objective 


The  objective  of  this  research  is  to  develop  a  computational  capability  for  large-scale  parallel 
simulations  of  microstructure  evolution  in  crystalline  solids.  The  task-oriented  goals  are  to 

(1)  design  a  flexible  computational  finite  element  (FE)  framework  for  phase-field  modeling; 

(2)  implement  this  framework  on  the  large-scale  computing  platforms  maintained  by  the  High 
Performance  Computing  Modernization  Office  (HPCMO);  and  (3)  perform  verification  and 
validation  (V&V)  simulations,  focusing  on  the  problem  of  mechanical  twinning  in  lightweight 
materials  of  interest  to  the  U.S.  Army  Research  Laboratory  (ARE).  Simulations  are  intended  to 
provide  new  insight  into  behavior  of  lightweight  metals  and  ceramics  concurrently  under 
development  in  the  Weapons  and  Materials  Research  Directorate  (WMRD)  for  vehicular  armor 
and  personnel  protection  applications. 


2.  Approach 


The  technical  approach  consists  of  the  following  two  steps:  (1)  development  of  a  new  phase-field 
theory  for  deformation  twinning  in  crystalline  solids  and  (2)  implementation  and  testing  of  the 
theory  in  EE  software.  Each  of  these  steps  is  summarized  in  a  separate  section  in  what  follows. 

2,1  Phase-field  Theory 

A  complete  description  of  the  phase-field  theory  developed  by  the  authors  will  be  available  in  a 
forthcoming  article  (1).  The  theory  is  summarized  below,  following  a  brief  introduction. 
Standard  notational  conventions  of  continuum  mechanics  are  used. 

Phase-field  models  have  been  successfully  applied  towards  a  variety  of  problems  in  materials 
science,  chemistry,  and  physics.  At  each  point  in  the  problem  domain,  one  or  more  phase-field 
variables  describe  the  state  of  the  substance.  The  term  “phase”  here  denotes  a  certain  crystal 
structure  or  lattice  configuration.  In  regions  of  uniform  phase,  order  parameters  acquire  discrete 
values  of  zero  or  one.  In  interfacial  regions,  order  parameters  enable  interpolation  between  pure 
phases.  In  addition  to  depending  on  usual  mechanical  and  thermal  state  variables  (e.g.,  strain 
and  temperature),  the  free  energy  density  of  a  substance  generally  depends  on  local  value(s)  of 
order  parameter(s)  and  spatial  gradients  of  order  parameter(s).  Surface  energies  of  phase 
boundaries  can  be  directly  associated  with  order  parameter  gradients  in  interfacial  regions. 

Cahn  and  Hilliard  (2)  and  Allen  and  Cahn  (5)  studied  thermodynamics  and  kinetics  of 
heterogeneous  material  systems  in  the  context  of  phase-field  models.  In  traditional  phase-field 
modeling,  it  is  typically  assumed  that  a  material  system  will  tend  to  evolve  towards  a  state  of 
minimum  free  energy,  subject  to  boundary  constraints  imposed  on  the  system.  Concepts  for 
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modeling  multiphase  systems  were  advaneed  by  Steinbaeh  et  al.  {4)  and  Steinbaeh  and  Apel  (5). 
Fried  and  Gurtin  (6)  and  Gurtin  (7)  developed  order  parameter  theories  ineorporating 
configurational  forces  in  the  context  of  geometrically  nonlinear  continuum  mechanics  and 
thermodynamics.  Review  articles  describe  numerical  techniques  and  applications  (5,  9). 

Mathematically,  a  scalar  order  parameter  field  is  denoted  by  ri(X,t) ,  with  X  material  coordinates 
and  t  time.  The  two  distinct  phases  represented  are  the  original  crystal  (i.e.,  the  parent)  and  the 
twin,  with  interfaces  between  phases  representing  twin  boundaries.  In  terms  of  the  order 
parameter, 

7  =  0  VX  e  parent;  7  =  1  VX  e  twin;  0  <  7  <  1  VX  e  interfacial  regions  .  (1) 

Let  X  =  x(X,t)  =  X  +  u(X,t)  denote  spatial  coordinates  of  a  material  particle,  with  u  the 
displacement.  The  deformation  gradient  F  is 

F  =  Vx,  (2) 

with  =  5/ the  material  gradient.  Deformation  gradient  (equation  2)  is  further 
decomposed  as 

F  =  F^F\  (3) 

where 

F®^(X,t)  =  FF'' '  :=  elastic  deformation;  F''[7(X,t)]  ;=  "stress-free"  twinning  shear.  (4) 

Specifically,  the  twinning  shear  is  interpolated  in  interfacial  regions  as  follows: 

F”  =  1  +  [^(7)]yoS  ®  m  ,  (5) 

where  (p  is  an  interpolation  function  obeying  0  <  ^  <  1,  ^(0)  =  0,  (p{l)  =  1,  ^  is  monotone.  A 
representative  function  also  satisfying  {dcp! drf)\^_^  =  {dcp! =  0  used  later  is  {10) 

^3  =  «7^ +2(2-«)7^ +(«-3)7'^,  (6) 

where  or  is  a  scalar  constant  within  the  limits  0  <  or  <  6 .  The  total  free  energy  functional  for  a 
body  of  reference  volume  i7  is  written  as 

4' =  \w{F,ri)dn  +  \f{ri,Vri)dn,  (7) 

n  n 

where  W  is  the  elastic  strain  energy  density  (generally  non-zero  in  parent,  twin,  and  interfacial 
regions),  and  where/ accounts  for  interfacial  energy  in  twin  boundary  regions.  Strain  energy  per 
unit  reference  volume  is  of  the  functional  form 

W  =  W[F^  (F,  7),  7] ,  E*"  =  (1  /  2)(C*'  -!)  =  (!/  2)(F*^^F''  - 1) ,  (8) 
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with  C®  the  elastic  deformation  tensor,  the  elastic  strain  tensor,  and  a  superposed  T  denoting 
the  transpose.  For  any  value  of  t] ,  strain  energy  density  vanishes  at  null  elastic  strain: 

W(0,-)  =  0  .  (9) 

A  quadratic  form  for  the  strain  energy  density  is  assumed.  Higher-order  elastic  coefficients  can 
be  incorporated  by  extension  (JJ).  Dependence  of  on  fj  manifests  explicitly  only  via 

anisotropic  elastic  coefficients.  Strain  energy  density  and  second-order  moduli  are  written 

1F  =  (1/2)E'^:C(;7):E^  C(77)  =  5'1F/(5E*^5E®)I  .  (10) 

IE  =0 

As  usual,  indices  =  ^badc  -  ^cdab  the  Laue  group  of  the  crystal  dictate  any  other 
symmetries  of  C(0)  in  the  crystallographic  frame  of  reference.  For  a  compound  twin  in  a 
centrosymmetric  structure,  reorientation  matrix  Q  associated  with  twinning  is  {12) 

Q  =  2m0m-1.  (11) 

Elastic  coefficients  of  the  fully  twinned  crystal  are  related  to  those  of  the  parent  by 

^ABCD  (1)  Qae  Qbf  Qcg  Qdh  ^efgh  (0).  (12) 

Elastic  coefficients  in  interfacial  regions  are  interpolated  the  same  way  as  the  twinning  shear: 

C(7)  =  C(0)  +  [C(1)-C(0)]^,  (13) 

where  (p{ri)  obeys  equation  6.  In  the  isotropic  elastic  approximation,  letting  A,  denote  Eame’s 
constant  and  //  the  shear  modulus,  the  elasticity  tensor 

^ABCD  ~  ^^AB^CD  +  K^aC^BD  +  ^AD^BC  )  ’  (  1 4) 

SO  that  C(0)  =  C(l)  and  W(E^  ,7])  -^W(E^)  does  not  explicitly  depend  on  7 .  The  local 
interfacial  energy  per  unit  reference  volume  follows  the  Cahn-Hilliard  formalism  (2): 

/(7,V7)  =  /o(7)  +  k:(V7  0V7),  (15) 

with  K  a  symmetric  second-order  tensor,  here  assumed  constant.  When  interfacial  energy  is 
isotropic,  k  =  /fl  and 

f{TJ,VTj)=f^{Tj)  +  K\VTj\\  (16) 

Prescribed  for  is  a  standard  “double-well”  potential  {6,  8): 

/o=V(l-7)%  (17) 

withH  a  non-negative  constant.  As  shown  in  reference  1,  in  the  isotropic  approximation  H  and 
K  are  related  to  equilibrium  energy  per  unit  area  F  and  thickness  /  of  an  unstressed  interface 
via 
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K  =  3ri/A,  A  =  \2r/i. 


(18) 


The  total  free  energy  funetional  ^  of  equation  7  beeomes,  using  equations  10  and  15: 

•F  =  (1  / 2) j E"  :  C(;7) :  E +  \[Ari\\-rif +K:{Vri®V ri)\dQ .  (19) 

n  n 

For  isotropic  elastic  and  interfacial  energies, 

'F  =  \[{XI2){ixE^f  ^  :^^]dQ  +  \[Ari\\-rif  +K\Vrif]dn .  (20) 

D  Q 

The  following  variational  equation  is  posited  that  will  suggest,  upon  application  of  Hamilton’s 
principle,  local  or  strong  forms  of  static  equilibrium  equations  and  boundary  conditions: 

dF-\  UdxdS  -  j  hdr/dS  =  0 ,  (21) 

dO  dQ 

where  t  is  a  mechanical  traction  vector  per  unit  reference  area,  dS  is  a  surface  element  of  Q ,  and 
/z  is  a  scalar  conjugate  force  to  variations  of  the  order  parameter.  Taking  the  first  variation  of  the 


interfacial  energy  and  applying  the  divergence  theorem, 

^  j  fdQ  =  j  (a/o  /  dri)dridV  -  J2k  :  [V{Vri)\dridn  +  j  2k  :  [{V  ri)®n\dridS ,  (22) 

O  D  Q  dQ 

with  n  the  unit  outward  normal  to  dQ .  Taking  the  first  variation  of  the  strain  energy  with  F  and 
7  independent, 

WdQ  =  j  (dJV/d7^)d7^dQ-j[V^(dJV/dF)]^dxdQ+  j  [n^(dJV / dF)]-dxdS  .  (23) 

Q  Q  Q  dQ 

It  follows  that  within  i7 ,  Euler-Lagrange  equations  are 

v.(aiE/aF)|^  =  v.p  =  o,  v,(aiE/av,vj|^=v,E^=o;  (24) 

a/o/a7-2K:[V(V7)]  +  (aif/a7)|,=0,  dfJdrj-2K,,V,V,rj  +  idW/dr])l=0-,  (25) 

where  P  is  the  first  Piola-Kirchhoff  stress  tensor.  Corresponding  boundary  conditions  are 

t  =  Pn,  (26) 

/z  =  2k  :  (V  7  0  n) ,  h  =  2/y^^n^V^7 .  (27) 

The  first  Piola-Kirchhoff  stress  also  obeys 

p  =  (aiF/aF)|^  =(aiF/aE*^):(aEVaF'^):(aFVaF)|^  =F''I:F''-^  (28) 
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The  elastic  second  Piola-Kirchhoff  stress  is 

'L  =  dW  IdE^  =  .  (29) 

The  partial  derivative  of  1F[E'^(F,;7),;7]  with  F  fixed  is  computed  as 

(dJr/dfj)i=(dW/dfj)L+(dW/dE^)i  :(dE^ldrj)\  .  (30) 

It  IF  Ip 

Furthermore, 

(dW!  dr])\^,  =  (1  /  2)F'' :  (5C  /  5^) :  F®  =  (1  /  2)(d(pl  d7])E^  :  [C(l)  -  C(0)] :  F’'  ,  (31) 

(dW/dEH  :{dE^ldri)\  =-{i::[C^{%®m)E''-^])y^{d(pl dri)  =  -Ty,{d(pl ^l^).  (32) 

\t}  if 

It  can  be  verified  that  r  is  a  kind  of  resolved  shear  stress  on  the  habit  plane  in  the  direction  of 
twinning  shear.  Phase  equilibrium  condition  (equation  25)  can  be  rewritten 

dfj dr] -2k\[W(W ri)'\  =  {(H 2)E^  \[C(9)-C(\)y.E^  +TY^}(d(pldri).  (33) 

In  the  isotropic  approximation,  choosing  cp  from  equation  6  and  from  equation  17, 
equation  33  becomes 

T[arj  +  '2(2-a)rj^  +2(a-'2)rj^]  =  (H Yy)[Arj(\-'2rj  +  2rj^)- .  (34) 


Both  sides  of  equation  34  vanish  identically  in  regions  of  uniform  phase  where  7  =  0  or  7  =  1. 

The  preceding  derivations  focus  on  the  geometrically  nonlinear  theory,  with  only  a  single  twin 
system  (i.e.,  a  single  order  parameter).  Geometric  linearization  of  the  theory,  extension  of  the 
theory  to  account  for  multiple  twin  systems  (i.e.,  multiple  order  parameters),  and  existence  of 
solutions  to  the  energy  minimization  problem  corresponding  to  equation  21  are  discussed  in 
reference  1. 

2.2  Finite  Element  Formulation 

The  FE  method  is  used  to  seek  solutions  for  equilibrium  or  local  minimum  energy  states  of  a 
body  subjected  to  boundary  conditions.  At  each  reference  point  X,  primary  solution  variables 
are  displacement  u  and  order  parameter  7 .  If  these  variables  and  requisite  material  properties  are 

known,  then  all  elastic  field  variables  and  interfacial  quantities  can  be  computed  via 
mathematical  operations  given  in  section  2.1. 

The  body  of  reference  volume  £2  is  discretized  into  a  number  of  standard  FEs  with  shape 
functions  (X) .  Let  u.(t)  and  rj^t)  denote  instantaneous  values  of  displacement  and  order 

parameter,  respectively,  at  node  i.  Displacement  and  order  parameter  fields  are  represented  in  a 
EE  context  as 


u(X,0  =  u,(0A,(X) ,  7(X,0  =  ri,(t)NyX)  , 


(35) 
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where  summation  proeeeds  over  all  nodes  z;  however,  only  those  nodes  supporting  the  element(s) 
containing  or  bounding  point  X  have  A^,  (X)  ^  0  .  Material  gradients  of  equation  35  follow  as 

F  =  l  +  Vu  =  l  +  u,.  0  VA^.  ,Vri  =  .  (36) 


The  FE  method  is  used  to  seek  solutions  of  global  or  weak  forms  of  the  equilibrium  conditions. 
Strong  (local)  forms  are  not  needed  by  the  numerical  algorithms.  Addressed  in  what  follows  are 
the  following  kinds  of  boundary  conditions: 

do  =  dOj^  n dOp', 

d^M=d^M,D^dn^^^,0  =  dn^^^ndn^y,  1 

>  mechanical  conditions  (37) 

u(X,t)  prescribed  on  dO^^  p,,  t(X,t)  =  0  on  dO^^  I 

dnp=dnppUdnpj„0  =  dnppndnpy,  1 

1  phase  held  conditions 

ri(X,t)  prescribed  on  dOp^,  h(X,t)  =  0  on  dOp  j^A 


Dirichlet  conditions  for  displacement  and  the  order  parameter  are  applied  on  ^  and  dOp  ^  , 
respectively.  Neumann  conditions  corresponding  to  free  surfaces  are  applied  on  ^  and 
d£2p  ^ ,  respectively.  The  present  treatment  can  be  extended  to  address  arbitrary  Neumann 

conditions  if  terms  accounting  for  external  work  are  incorporated  in  the  total  energy  functional 
whose  stationary  points  are  sought. 


The  free  energy  functional  is  equation  19  in  the  general  anisotropic  case.  Nodal  equilibrium 
conditions  are  obtained  by  substituting  equations  35  and  36  into  equation  19  and  differentiating 
with  respect  to  nodal  degrees  of  freedom: 

0  =  a  'A/au.  =  j  FVN.dn  =  j  [F*'(C  :  E'')F’'-^]VA,Jf2 ,  (38) 

n  Q 


0  =  a  “A  / az/,.  =  j  {2.fz7(l  -  3;/  +  2;/")}  N^dO  +  j  {2kV 

Q  Q 

+  j{(l/2)E*^:[C(l)-C(0)]:E*^  -TYQ^{d(p  !  dri)Njd£2. 

n 


(39) 


A  conjugate  gradient  algorithm  is  used  to  seek  minimum  energy  states  corresponding  to 
equations  38  and  39,  following  prescription  of  initial  conditions. 

2,3  Phase-field  Software  Framework 

The  above  FE  numerical  formulation  is  implemented  into  a  standalone  software  framework.  In 
contrast  with  more  commonly  employed  approaches  relying  on  a  single  executable  to  handle  all 
variability  in  problems  and  data,  our  computational  framework  enables  users  to  rapidly  develop 
applications  tailored  to  their  particular  needs.  This  flexibility  is  achieved  through  modern 
software  development  design  and  practices.  The  core  computational  kernels  (e.g.,  FE 
methodology,  quadrature  integration  rules,  phase-field  energy  densities,  and  minimization 
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solvers)  are  implemented  in  ANSI  C++,  which  affords  great  code  portability,  modularity,  and 
code  reuse.  In  contrast,  high-level  user-interfaces  are  accessible  both  from  C++  and  Python  for 
rapid  application  development.  The  modular  structure  of  the  framework  greatly  simplifies 
extension  efforts,  for  example,  adding  new  models  of  material  or  interface  behavior.  Moreover, 
such  structure  allows  for  a  wide  spread  use  of  generic  programming  concepts  that  enable 
separation  of  algorithms  from  the  data  on  which  they  operate.  By  way  of  example,  we 
implemented  a  family  of  generic  nonlinear  minimization  solvers  (Jacobi-preconditioned 
conjugate-gradient  and  Newton-Raphson)  that  can  be  employed  to  solve  problems  derived  from 
the  FE  method,  as  well  as  those  characterized  by  an  explicit  function.  Input  and  output  are 
handled  exclusively  via  ARL’s  Extensible  Data  Model  and  Eormat  (XDME)  (25).  Input  consists 
of  a  EE  domain  discretization  (mesh)  along  with  model  parameters  contained  in  an  XDME  file. 
Owing  to  the  fact  that  XDME  contains  filters  from  popular  computer-aided  design  (CAD)  mesh 
formats,  external  off-the-shelf  meshing  tools  can  be  employed.  An  output  file  contains  a 
snapshot  of  the  entire  internal  state  of  a  phase-field  application  and  thus  can  be  used  to  resume 
its  execution.  In  addition,  the  output  fide  may  be  used  for  visualization  purposes  as  many 
visualization  packages  are  able  to  directly  read  XMDE  files. 

The  phase-field  software  framework  can  take  advantage  of  large-scale  distributed  memory  high- 
performance  computers  (HPCs)  by  recourse  to  the  Message  Passing  Interface  (MPI)  (26).  Our 
implementation  relies  on  conventional  non-overlapping  domain  decomposition  methodology 
(27)  in  which  a  PE  mesh  is  distributed  among  processors  and  compatibility  conditions  are 
enforced  across  sub-domain  boundaries.  Mesh  decomposition  is  carried  out  via  a  set  of  XDME 
tools  that  perform  mesh  partitioning  as  well  as  establish  mappings  between  shared  mesh  vertices 
on  all  sub-domain  boundaries.  The  main  advantage  of  the  non-overlapping  domain 
decomposition  methodology  is  it  relative  simplicity  and  limited  scope  of  modifications  to  serial 
implementations  while  still  offering  remarkable  performance  and  scalability  on  modern  HPC 
systems. 


3.  Results 


3.1  Twin  Nucleation  in  Magnesium 

Homogeneous  twin  nucleation  refers  to  creation  of  a  twin  embryo  within  an  otherwise  perfect 
single  crystal  (72).  Analytical  models  based  on  free  energy  variations  in  the  context  of  phase 
transformations  have  been  applied  to  describe  twin  nucleation  (72,  13).  Such  approaches 
consider  nucleation  of  a  twin  embryo  of  idealized  geometry — an  elliptical  cylinder  in  two 
dimensions  or  an  ellipsoid  in  three  dimensions — embedded  in  an  infinite  medium,  with  a 
perfectly  bonded  sharp  interface  separating  inclusion  from  surrounding  medium  (i.e.,  the  matrix). 
The  solution  technique  involves  simultaneous  solution  of  two  equilibrium  equations  associated 
with  stationary  points  of  the  Gibbs  free  energy  change  associated  with  twinning.  These  two 
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equations  yield  the  eritieal  size  and  aspeet  ratio  of  the  inelusion,  as  outlined  in  referenees  1  and 
13.  Apparently,  exaet  solutions  are  available  only  for  the  approximations  of  linear  elastieity, 
isotropie  surfaee  energy  assoeiated  with  the  twin  boundary,  and  traetion  boundary  eonditions  at 
infinity.  The  solution  provides  the  eritieal  size  and  shape  of  an  inelusion  for  a  given  set  of 
material  properties — elastie  eonstants,  twin  boundary  surfaee  energy,  and  twinning 
transformation  shear — and  far-field  applied  stress.  At  the  eritieal  aspeet  ratio,  the  eritieal  size 
eorresponds  to  unstable  equilibrium,  i.e.,  a  saddle  point  on  the  free  energy  surfaee.  At  a  given 
far-field  applied  stress,  a  twin  nueleus  smaller  than  the  eritieal  size  will  tend  to  shrink  and 
disappear,  while  one  larger  than  the  eritieal  size  will  tend  to  grow  in  an  unstable  manner.  The 
larger  the  applied  stress  eomponent  resolved  on  the  habit  plane  in  the  twinning  direetion,  the 
smaller  the  eritieal  size,  meaning  that  nueleation  of  a  twin  of  fixed  size  beeomes  more 
energetieally  favorable  as  the  resolved  shear  stress  inereases. 

The  phase-field  approaeh  is  used  in  the  present  work  to  model  the  problem  of  homogeneous  twin 
nueleation.  Numerieal  results  obtained  from  the  phase-field  model,  under  assumptions  of 
geometrie  linearity  and  plane  strain  boundary  eonditions,  are  eompared  with  the  analytieal 
solution  {13).  This  exereise  provides  validation  for  the  phase-field  model  and  numerieal 
methods  advaneed  here.  Additional  numerieal  simulations  eonsider  effeets  of  anisotropie  twin 
boundary  surfaee  energy,  variable  equilibrium  thiekness  of  the  twin  boundary  region,  various 
habit  plane  direetions  (i.e.,  lattiee  orientations),  and  various  boundary  eonditions  (Diriehlet  vs. 
Neumann,  shear  vs.  tension,  and  domain  shapes).  These  additional  faetors  eannot  all  be 
addressed  via  known  analytieal  elastieity  solution  teehniques. 

Pure  magnesium  (Mg)  single  erystals  are  studied.  Mg,  a  metallie  erystalline  solid  at  standard 
temperature  and  pressure,  is  of  partieular  interest  for  lightweight  vehieular  proteetion 
applieations  beeause  of  its  low  mass  density.  Properties  are  listed  in  table  1  with  supporting 
referenees.  Mg  exhibits  a  hexagonal  erystal  strueture  and  is  eentrosymmetrie,  with  e/a  ratio 
1.6235.  The  twinning  system  of  eonsideration  is  the  primary  one;  <  10  1 1  >  {1012}  with 
twinning  shear  =  (3  -c^  /  a^)/(3'^^c/a)  =  0.1295  .  Voigt  averages  are  used  for  isotropie  elastie 

eonstants.  Mg  single  erystals  are  nearly  isotropie  elastieally  [Ch^Cbs,  Ci2«Ci3,  and  C44  «(Cii  - 
Ci2)/2],  so  this  isotropie  elastie  approximation  is  reasonable. 


Table  1.  Properties  for  pure  Mg  single  crystals. 


Parameter 

Value 

Definition 

Reference 

To 

0.1295 

Shear  for  <  10  1 1  >  {1012}  twin 

12 

T 

24.0  GPa 

Lame  constant  (Voigt  average,  0  K) 

14 

19.4  GPa 

Shear  modulus  (Voigt  average,  0  K) 

14 

r 

117  mJ/m^ 

Twin  boundary  energy  (0  K) 

15 

i 

1.0  nm 

Equilibrium  boundary  thickness 

16 
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Sought  are  critical  applied  strains  for  twin  nucleation  for  the  simulation  cases  listed  in  table  2, 
which  encompass  various  boundary  conditions,  material  properties,  and  lattice  orientations. 
Here,  L/H  denotes  the  aspect  ratio  of  the  simulation  domain,  where  //is  fixed  at  50  nm.  Critical 
strains  are  obtained  in  practice  by  seeding  the  domain  with  an  inclusion  of  initial  radius 
Oq  =  3  nm  (which  is  on  the  order  of  the  minimum  stable  size  observed  from  first  principles 

calculations  [/5]),  as  shown  in  figure  1,  and  then  applying  the  boundary  deformation  in  small 
increments  of  magnitude  0.001.  For  each  increment,  system  degrees  of  freedom  relax  in 
conjunction  with  energy  minimization,  performed  numerically  via  the  conjugate  gradient 
technique.  The  minimum  applied  strain  for  which  decay  in  size  of  the  twin/inclusion  does  not 
occur  is  deemed  the  critical  strain  .  The  corresponding  average  shear  stress  is  . 


Table  2.  Phase-field  simulation  parameters,  predicted  critical  strain  for  twin  nucleation  y^,  and  average 
equilibrium  shear  stress  after  twin  has  fully  grown. 


Case 

0rad 

1  nm 

L/H 

Boundary 

Loading 

rc 

Tc/H 

1 

0 

1 

1 

1 

Dirichlet 

shear 

0.050 

0.020 

2 

0 

1 

1 

1 

Neummann 

shear 

0.050 

0.000 

3 

0 

1 

1 

1.5 

Dirichlet 

shear 

0.049 

0.012 

4 

0 

1 

4“ 

1 

Dirichlet 

shear 

0.046 

0.023 

5 

0 

0.1 

1 

1 

Dirichlet 

shear 

0.049 

0.044 

6 

71/6 

1 

1 

1 

Dirichlet 

shear 

0.050 

0.042 

7 

71/4 

1 

1 

1 

Neummann 

tension 

0.045 

0.005 

“Anisotropic  surface  energy:  kJ2  =  =  k. 


Figure  1 .  FE  mesh  and  initial  conditions 
of  an  undeformed  shape. 

Figure  2  shows  equilibrium  order  parameter  contours  for  each  simulation  at  the  applied  critical 
strain.  Figure  3  shows  corresponding  shear  stresses.  For  cases  1-6,  the  local  shear  stress  shown 
is  (Jy2 ;  for  case  7,  the  stress  shown  is  that  acting  on  a  45°  plane,  {a 22  -crjj)/2 .  Results  are 

mesh- insensitive  except  for  those  of  case  5,  wherein  growth  of  the  twin  is  confined  to  the  center 
of  the  domain  where  the  mesh  is  refined  enough  to  resolve  gradients  of  the  order  parameter  at  the 
interface.  Predictions  of  critical  nucleation  strain  from  the  phase-field  simulations  are  compared 
with  predictions  of  other  models  in  table  3. 
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Figure  2.  Equilibrium  order  parameter  field  (fully  formed  twin)  at  critical  strain  for  simulations  listed  in  table  2: 
(a)  case  1,  (b)  case  2,  (c)  case  3,  (d)  case  4,  (e)  case  5,  (f)  case  6,  and  (g)  case  7. 


0.87  _ 

0.55  n 

SiH 

0.43  ■ 


(C) 


(d) 


(e) 


(f) 


(g) 


Figure  3.  Equilibrium  effective  shear  stress  field  (fully  formed  twin)  at  critical  strain  for  simulations  listed  in  table  2: 
(a)  case  1,  (b)  case  2,  (c)  case  3,  (d)  case  4,  (e)  case  5,  (f)  case  6,  and  (g)  case  7. 

Table  3.  Critical  strain  for  twin  nucleation  or  growth:  comparison  with  other  models. 


Model 

Reference 

Remarks 

yc 

Phase-field 

This  work 

Cases  1-7,  Table  2 

0.045-0.050 

Elastic  inclusion 

13 

Analytical  solution,  infinite  medium 

0.038 

Theoretical  strength 

17 

=  YoiiiTt) 

0.021 

Atomic  calculation 

18 

Glide  of  1/17[1  oil] (To  12)  partial 

0.004 

10 


3,2  Twin  Nucleation  in  Calcite 

Considered  next  is  the  problem  of  twin  nucleation  in  calcite  (CaCOa).  CaCOs  is  of  interest  for  a 
number  of  reasons.  Experiments  in  the  1930s  through  the  1960s,  conducted  primarily  in  the 
Soviet  Union  {19,  20),  demonstrated  the  phenomenon  of  elastic  twinning  in  CaCOa  single 
crystals  subjected  to  concentrated  surface  loading,  e.g.,  indentation.  Elastic  twins  partially  or 
fully  disappear  upon  load  removal.  This  phenomenon  is  realistically  modeled  by  the  present 
phase  field  approach  that  considers  twinning  in  the  non-dissipative  and  null  temperature  (0  K) 
limit.  CaCOa  is  an  ideal  material  for  the  study  of  twinning  because  it  is  transparent  (hence  twins 
are  easily  observed)  and  exhibits  little  plasticity  (i.e.,  uncorrelated  slip  of  dislocations)  or 
fracture,  relative  to  twinning,  at  low  temperatures.  Requisite  properties  of  CaCOa  are  listed  in 
table  4  with  supporting  references. 

Table  4.  Properties  for  pure  CaCOs  single  crystals. 


Parameter 

Value 

Definition 

Reference 

ro 

0.694 

Shear  for  twin 

21 

A 

61.7  GPa 

Lame  constant  (Voigt  average,  0  K) 

22 

40.2  GPa 

Shear  modulus  (Voigt  average,  0  K) 

22 

r 

183  mJ/m^ 

Twin  boundary  energy  (0  K) 

23 

i 

1.0  nm 

Equilibrium  boundary  thickness 

23 

The  problem  of  study  is  illustrated  schematically  in  figure  4(a).  A  single  crystal  oriented  for 
twinning  is  subjected  to  concentrated  surface  loading  via  a  wedge-shaped  indenter.  The  indenter 
is  assumed  to  adhere  rigidly  to  the  crystal  surface.  The  crystallography  and  coordinate  system 
are  discussed  in  more  detail  in  reference  1 1 . 


Figure  4.  Indentation  of  CaC03  single  crystal  oriented  for  twinning:  (a)  sketch  of  boundary  value  problem  {19) 
and  predicted  order  parameter  field  from  phase-field  simulations:  (b)  90°  wedge  and  (c)  120°  wedge. 


Contours  of  the  order  parameter  field  predicted  from  phase-field  simulations  are  shown  in 
figures  4(b)  and  4(c).  The  indentation  depth  is  4.75  nm.  Notice  that  in  each  case  the  major  twin 
forms  asymmetrically  under  only  one  side  of  the  indenter  because  the  shear  deformation  on  the 
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opposite  side  is  directed  in  the  anti-twinning  sense.  The  twin  shapes  are  long  and  thin,  with 
sharp  tips,  as  observed  in  experiments  {19,  20).  It  is  found  that  the  threshold  for  twin  nucleation 
is  highly  sensitive  to  the  presence  of  initial  defects  in  CaCOa,  and  that  the  length  L  of  the  twin 
increases  monotonically  with  increasing  load  P.  Both  such  trends  in  results  have  also  been 
observed  experimentally  {19,  20,  24).  Also  noteworthy  are  the  small  secondary  twins  that  appear 
spontaneously  on  the  right  side  of  the  indenter,  at  the  free  surface. 

Results  shown  in  figures  4(b)  and  4(c)  are  obtained  using  the  phase  field  theory  under  the 
approximation  of  linear  isotropic  elasticity.  A  more  complete  study  of  twinning  in  CaCOa  single 
crystals  using  the  phase  field  approach  will  appear  in  a  forthcoming  publication  {11),  which  will 
also  incorporate  nonlinear  and  anisotropic  elasticity.  CaCOa  is  of  Army  interest  because  it  is  a 
lightweight  nonmetallic  mineral  with  ductility  (i.e.,  a  kind  of  “soft  ceramic”)  whose  study  may 
provide  qualitative  insight  into  inelastic  phenomena  that  arise  in  other  transparent  ceramics  (used 
for  armor  applications)  that  undergo  twinning  such  as  alumina  (AI2O3)  and  aluminum  oxynitride 
(AlON).  CaCOa  belongs  to  the  same  crystallographic  space  group  as  a-alumina  (i.e.,  sapphire); 
both  belong  to  centrosymmetric  point  groups  with  trigonal  symmetry. 

3,3  Parallel  Performance 

We  investigated  parallel  performance  of  the  phase-field  software  framework  for  the  problem  of 
twin  nucleation  in  magnesium  (see  figure  1).  Our  measurements  were  performed  on  MJM  HPC 
system  located  at  ARL’s  Department  of  Defense  (DoD)  Supercomputing  Resource  Center. 

The  first  performance  test  was  focused  on  the  measurement  of  parallel  speedup.  Parallel 
speedup,  is  commonly  defined  as 
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p 


(40) 


where:  7^ ,  denote  the  execution  time  of  the  sequential  and  parallel  algorithm  on  p  processors, 
respectively.  We  measured  for  two  discretizations  of  the  computational  domain  containing 

257,094  and  502,780  linear  triangular  elements  and  the  number  of  processors  in  the  1  to  2048 
range.  The  speedup  as  a  function  of  the  numbers  of  processors  is  plotted  in  figure  5(a)  along 
with  the  ideal  speedup  {S^=  p).  Close  to  ideal  parallel  speedup  may  be  observed  in  the  1  to  512 

processor  range.  Outside  this  range  (for  processor  counts  larger  than  512)  the  speedup  is 
markedly  reduced.  This  behavior  is  expected  as  the  average  number  of  elements  per  processor 
for  1024  processors  is  reduced  to  merely  251  and  490,  respectively.  Under  such  circumstances 
the  computational  expense  becomes  significantly  influenced  by  the  cost  of  inter-processor 
communication. 
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Figure  5.  (a)  Measured  values  of  parallel  speedup  vs.  number  of  processors  and  (b)  execution  time  per  one 
conjugate  gradient  iteration  vs.  number  of  processors. 


Apart  from  the  parallel  speedup  measurement,  we  also  measured  the  execution  time  under 
condition  of  constant  number  of  elements  per  processor.  To  this  end,  we  initially  created  a  mesh 
containing  6,192  linear  triangular  elements.  Subsequently,  this  initial  mesh  was  uniformly 
subdivided  five  times  yielding  a  sequence  of  meshes  containing  24,768;  99,072;  396,288; 
1,585,152;  and  6,340,680  triangular  elements.  It  bears  emphasis  that  each  uniform  subdivision 
step  raises  the  number  of  elements  by  4.  Thus,  when  the  number  of  processors  is  increased 
accordingly  by  4  with  each  subdivision  the  number  of  elements  per  processor  remains,  on 
average,  unchanged.  We  measured  the  execution  time  for  3,000  iterations  of  Jacobi- 
preconditioned  nonlinear  conjugate  gradient.  The  execution  time  per  iteration  is  plotted  in 
figure  5(b).  As  expected,  the  execution  time  is  the  smallest  for  the  sequential  run  and  increases 
significantly  for  parallel  runs.  We  emphasize,  however,  that  the  parallel  execution  times  are 
within  18%  of  each  other,  a  clear  indication  of  very  good  parallel  performance. 


4.  Conclusions 


Phase-field  theory  and  simulation  software  have  been  developed.  The  theory  accounts  for 
deformation  twinning  in  deformable  crystalline  solids,  with  equilibrium  equations  obtained  via  a 
variational  principle  in  the  null  temperature  limit.  Numerical  solutions  to  weak  forms  of 
governing  equations  are  obtained  via  a  conjugate  gradient  solver  and  the  FE  method. 

Two  fundamental  problems  in  materials  physics  have  been  investigated.  The  first  involves 
homogeneous  nucleation  of  a  twin  in  a  Mg  single  crystal.  Critical  shear  strains  for  nucleation 
obtained  numerically  using  the  phase-field  approach  are  in  reasonable  agreement  with  those 
obtained  analytically  in  the  sharp-interface  limit.  The  second  involves  twin  nucleation  in  a 
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CaCOs  single  crystal  subjected  to  indentation  loading.  Long,  thin,  asymmetric  twins  with  sharp 
cusp-like  tips  have  been  observed  in  numerical  simulations,  and  are  in  qualitative  agreement  with 
experimental  observations. 

All  results  obtained  are  predictive.  The  theory  and  numerical  implementation  rely  only 
fundamental  material  properties  (e.g.,  elastic  constants  and  twin  boundary  surface  energies)  that 
can  be  obtained  independently  through  physical  experiments  and/or  first-principles  calculations. 
Use  of  the  model  and  software  does  not  require  any  ad-hoc  parameter  fitting. 
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6.  Transitions 


6.1  Journal  Articles 

Clayton,  J.  D.;  Knap,  J.  A  Phase  Field  Model  of  Deformation  Twinning:  Nonlinear  Theory  and 
Numerieal  Simulations.  Aeeepted  for  publieation  in  Physica  0,2011. 

Clayton,  J.  D.;  Knap,  J.  Phase  Field  Modeling  of  Indentation- indueed  Twinning  in  Transparent 
Crystals.  I.  Linear  Theory,  in  preparation,  2011. 

Clayton,  J.  D.;  Knap,  J.  Phase  Field  Modeling  of  Indentation- indueed  Twinning  in  Transparent 
Crystals.  II.  Nonlinear  Theory,  in  preparation,  2011. 

6.2  Invited  Presentations 

Clayton,  J.  D.;  Knap,  J.  Phase  Field  Modeling  of  Twin  Nucleation  in  Magnesium  and  Calcite 
Single  Crystals,  Invited  presentation.  Center  for  Advaneed  Vehicular  Systems,  Mississippi 
State  University,  Starkville,  MS,  September  28,  2010. 

6.3  Transition  to  ARL  Mission  Programs 

This  work  will  tentatively  continue  to  receive  support  under  the  fiscal  year  2011  (FY 11)  WMRD 
Innovation  Project  program  under  the  project  titled  “Twinning  in  Transparent  Materials,”  which 
also  involves  M.  Greenfield,  J.  Swab,  and  C.  Hilton  of  WMRD.  Analytical  modeling,  phase- 
field  simulations,  and  physical  experiments  are  planned  that  will  probe  twinning  phenomena  in 
calcite  and  alumina  single  crystals. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


A1203 

alumina 

AlON 

aluminum  oxynitride 

ARL 

U.S.  Army  Research  Faboratory 

CaCOs 

calcite 

CAD 

computer-aided  design 

DoD 

Department  of  Defense 

FE 

finite  element 

FYll 

fiscal  year  2011 

HPCMO 

High  Performance  Computing  Modernization  Office 

HPCs 

high-performance  computers 

Mg 

magnesium 

MPI 

Message  Passing  Interface 

V&V 

verification  and  validation 

WMRD 

Weapons  and  Materials  Research  Directorate 

XDMF 

Extensible  Data  Model  and  Format 
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